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TECHNICAL  PROGRESS  REPORT 


Contract  Number  N00014-85-C-0852 


Radar  Cross-Section  of  Damped  Cylinders 
and  Dielectric  Strips 


Most  of  the  basic  logic  in  our  TDFD  RCS  code  was  in  place  at  the  begin¬ 
ning  of  this  quarter.  Consequently,  we  have  now  directed  our  primary 
concern  to 


1)  Improving  the  computer  efficiency  of  the  code. 

2)  Making  changes  which  lead  to  increased  accuracy. 

3)  Examining  a  second  canonical  problem,  a  thin  dielectric  strip,  to 
verify  further  the  code  accuracy. 

4)  Making  the  code  more  user-friendly. 

5)  Documenting  the  TDFD  RCS  code. 


By  way  of  Item  2  above,  we  have  recomputed  the  RCS  of  a  perfectly 

conducting  rod  . 5  m  in  radius,  bare  and  covered  by  a  damper  .5  m  thick.  The 

.  3 

damper  is  characterized  as  e/e0  -  /i/^o  -  1 ,  a  -  eo*/p  -  4  x  10  .  These 
values  were  selected  to  give  a  skin  depth  on  the  order  of  damper  thickness 
at  the  frequencies  (50  -  500  MHz)  for  which  the  calculations  were  run.  The 
TDFD  code  utilized  square  cells  4  cm  on  a  side  or  25  cells  to  a  cylinder 
diameter. 


While  our  TDFD  code  only  treats  the  TM  case,  we  are  able  to  simulate  TE 
problems  by  use  of  duality.  In  particular,  the  above  -  described  case  was 
rerun  with  the  cylinder  perfectly  magnetically  conducting  (a*  -  a  =  0) , 
and  the  damper  unchanged.  The  RCS  obtained  in  this  manner  for  TM  illumina¬ 
tion  is  the  same  as  one  would  obtain  for  TE  illumination  of  an  ordinary 
electrically  conducting  rod. 

Figures  1  and  2  show  our  TDFD  RCS  calculations  for  the  damped  and 
undamped  electrically  conducting  rod.  Overlaid  on  these  figures  are 
analytic  solutions  to  the  same  problem  as  determined  by  expansion  in 
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c  Rod  Radar  Cross  Section 


cylindrical  harmonics  (see  the  previous  quarterly  report  for  details  of  the 
harmonic  expansion).  These  results  should  be  compared  with  Figures  1  and  3 
of  the  previous  quarterly  report.  As  one  can  see,  the  frequency - doma i n 
ripple  which  had  been  present  on  the  response  of  the  damped  rod  is  now  gone. 
This  ripple  had  occurred  because  the  outer  boundary  of  the  problem  space  was 
somewhat  reflective.  Hence,  the  true  solution  was  "beating"  with  a 
coherent,  but  specious,  echo.  We  have  now  found  a  way  to  implement  a  much 
less  reflective  outer  boundary. 

(Unlike  3D  problems,  we  have  empirically  learned  that  2D  outer  bound¬ 
aries  are  least  reflective  if  a  pure  radiating  condition, 

f(r.  t)  -►  g(t  -  r/c)/7r 

is  used  without  a  damper  inside  the  boundary.  In  conjunction  with  this 
outer  boundary,  the  2D  code  should  not  be  run  with  At  at  the  Courant  limit, 
As/(/2c) ,  but  at  As/(2c) .  This  empirically  discovered  combination  gives 
less  outer-boundary  reflection  than  an  impedance  boundary,  a  damped  bound- 
arv,  a  pure  r  boundary,  a  Mur  boundary,  or  any  other  mixture  of  the 
above  .  ) 

Figures  3  and  4  show  our  TDFD  RCS  calculation  for  the  TE  simulation. 
Again,  these  results  are  overlaid  with  the  cylindrical-harmonic  solution. 
These  curves  correspond  to  Figures  2  and  4  of  the  previous  quarterly  report. 
Again,  the  new  results  may  be  seen  to  be  much  "cleaner". 

This  improvement  was  obtained  using  the  same  techniques  as  the  TM  case, 
but  with  an  additional  twist.  In  retrospect,  it  is  obvious  that  a  perfectly 
conducting  infinite  rod  illuminated  by  a  unipolar  EMP  polarized  along  the 
rod  axis  will  carry  a  d.c.  response  forever.  If  the  EMP  calculation  of  this 
response  is  stopped  suddenly  at  t0 ,  the  frequency  transform  of  the  result 
will  be  the  true  transform  f requency-convolved  with  sin(wt0/2)/(u>t0/2) .  We 
are  now  running  the  TE-like  TDFD  calculation  long  enough  that  2/t0  is 


axis 


smaller  than  any  angular  frequency  of  interest  and  then  turning  the  code  off 
"softly",  say  over  425  cycles. 

Additionally,  we  now  have  results  for  scattering  from  dielectric 
cylinders  .5  m  in  radius  and  characterized  by  t/c0  “  2  or  9.  (There  is  a 
fundamental  difference  between  these  two  cases.  At  e  ■=  2e0,  the  optically 
shortest  path  from  one  side  of  the  cylinder  to  the  other  is  through  the 
diameter.  At  £  **  9e0,  it  is  around  the  circumference.)  Figures  5-7  show 
overlays  of  these  calculations  and  the  harmonic  expansion.  It  may  be  seen 
that  the  overlays  begin  to  diverge  above  300  MHz,  especially  for  the  e  —  9 
case.  This  corresponds  to  about  8  cells  per  wavelength. 

Consequently,  the  dielectric  cylinder  calculations  were  rerun  with  As 
reduced  from  4  cm  to  2  cm.  Figures  8-10  show  the  overlays  based  on  the 
reduced  As.  These  new  results  are  considerably  improved,  especially  at  the 
higher  frequencies. 

It  is  our  belief  that,  although  great  improvement  has  been  achieved 
over  the  past  quarter,  we  are  still  quite  far  from  achieving  the  ultimate 
accuracy  TDFD  can  provide.  We  especially  hold  this  opinion  with  respect  to 
eliminating  the  broad  frequency-domain  ripple  appearing  on  our  latest  cal¬ 
culations  for  the  bare-rod  TE  simulation. 

Figures  11  and  12  present  overlays  for  the  second  canonical  problem, 
Item  3,  TM  scattering  off  a  dielectric  strip  .025  m  thick,  2  m  wide,  and 
characterized  by  a  «■  0 ,  £  —  2.  This  data  represents  the  monostatic  RCS  for 
a  45°  angle  of  illumination  with  respect  to  the  strip's  major  faces.  The 
TDFD  curves  in  these  figures  are  based  on  4  cm  square  cells.  At  the  higher 
frequencies,  2  cm  cells  would  doubtless  improve  the  code-code  agreement. 
The  frequency-domain  calculation  is  not  strictly  canonical,  but  is  actually 
a  modified  Galerkin  solution  with  three  complex  waves  included  to  represent 
the  strip  behavior  in  the  width  direction.  This  frequency-domain  solution 
is  due  to  Richmond,  and  was  brought  to  our  attention  by  Dr.  Stovall. 

With  respect  to  Item  4,  we  have  now  implemented  some  user-friendly- 
graphics  so  the  worker  has  a  visual  representation  of  the  problem  being 
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Dielectric  Rod  Radar  Cross  Section 


Echo  width  of  a  dielectric  rod  with  E  transverse 
to  the  rod  axis.  Based  on  rod  diameter  of  1  m 


Dielectric  Rod  Radar  Cross  Section 


Echo  width  of  a  dielectric  rod  with 
E  transverse  to  the  rod  axis.  Based 
on  a  rod  diameter  of  1  m  and  cells 


solved.  For  example,  Figure  13  illustrates  the  cell  occupancy  by  conductor 
or  damper  for  the  damped  electrically  conducting  rod  upon  which  Figures  1 
and  2  are  based. 

Lastly,  there  is  Item  5,  the  matter  of  TDFD  RCS  code  documentation. 
The  remaining  pages  of  this  report  constitute  our  first  draft  of  this  task. 


Architecture  of  the  Time-Domain  Finite-Difference  RCS  Code 


Equations  to  be  Solved 


We  assume  in  this  study  that  electromagnetically  linear  conditions 
prevail.  Then  the  total  field  can  be  separated  into  an  incident  field 
(which  would  be  the  field  in  the  absence  of  the  scatterer)  and  a  scattered 
field  (the  field  modification  caused  by  the  scatterer 's  presence): 


(ET,  HT)  =  (EinC,  ^ 


HinC) 


.  .scat 

(£  .  ^ 


Hscat) 


(1) 


Under  the  linearity  assumption,  both  the  incident  and  the  scattered  fields 
individually  satisfy  Maxwell's  equations. 

If  some  background  dissipation  (o^,  a^)  is  present,  the  incident  fields 
will  obey 
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In  the  presence  of  an  anisotropic  scatterer  with  frequency  dependent 
properties,  the  total  fields  conform  to 
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Subtraction  of  the  first  pair  of  equations  from  the  second  leaves  us  with 
the  version  of  Maxwell's  equations  obeyed  by  the  scattered  fields: 
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It  is  convenient  to  represent  the  inhomogeneous  parts  (or  incident 
parts)  of  eqs .  (6)  and  (7)  as 


[J 


T* 


1°*] [HS]  +  [J*]  +  [V  x  ESCat] 


[A  ]  =  [o0 ]  "(K  ] 


M  =  [^o]  1[«00  -  fo] 
[R]  -  [o0]'1[°r  -  ah] 

[A]  =  Ia0]'V] 


Method  of  Solution 


For  now,  we  shall  only  concern  ourselves  with  the  2D  TM  case.  Thus, 

E  ,  E  and  H  will  be  present.  Equation  (18)  reduces  to  a  scalar  equation, 
x  y  z 

but  (19)  remains  a  two  -  component  vector  equation. 

Figure  14  illustrates  a  typical  TDFD  2D  unit  cell,  and  shows  where  the 
three  field  components  associated  with  that  cell  are  located. 

In  Appendix  1,  we  will  derive  a  technique  for  solving  this  equation 
system  using  first-order  exponential  differencing.  The  result  of  this 
appendix  is  that,  if  we  omit  frequency  dependence,  HZ(I,J)  is  advanced 
according  to 


HZ ( I , J )  1  x  -  e  HZ(I ,  J)  -  (1  -  e 


-1  **- 
-u  a0At  . 


vn+l  _  0  tnH  ^co  o  -ljT~(][  j^n+1/2  (20) 


and  [E(I,J)]  is  advanced  according  to 


n+1/2  "  ^  ecJ  [£7o]At  n-1/2 

[E(I,J)]n+i/^  -  e  [E(I,J)]n  L/* 
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(21) 


(The  complication  of  frequency  dependence  will  be  considered  later.) 


As  we  have  said,  [ct0]  is  permitted  to  be  anisotropic.  In  the  actual 
code,  it  is  represented  by  a  total  of  five  arrays: 
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,  1  -  «  *  1  ....  -  ' 


SGX(I,J)  is  the  bulk  value  characterizing  the  xx  component  of  the 

conductivity  tensor  at  cell  (I,J) 

SGY(I,J)  is  the  bulk  value  characterizing  the  yy  component  of  the 

conductivity  tensor  at  cell  (I,J) 

SGXY(I.J)  is  the  hulk  value  characterizing  the  xy  and  yx  components  of 
the  conductivity  tensor  at  cell  (I,J)  (gyrotropic  materials  are  not 
permitted  in  the  present  code) 

SGCX(I,J)  is  the  surface  conductivity  in  the  x  direction  on  the  y- 

facing  surface  of  cell  (I,J) 

SGC Y(I,J)  is  the  surface  conductivity  in  the  y  direction  of  the  x- 

facing  surface  of  cell  (I,J) 

Thus,  the  actual  conductivity  seen  by  at  cell  (I,J)  is  given  by 

cr0(I,J)  ■=  (SGX(I.J-l)  +  SGX(I,J))/2  +  SGCX(I.J)  (22) 


a0(I,J)xy  =  (SGXY(I ,  J  - 1)  +  SGXY’CI ,  J)  )/2 
C70(I,J)  -  (SGXY  (I.J-l)"1  +  SGXY(  I ,  J ) " 1 )  ’ 1  •  2 

^o(I.J)  -  (  SGY  ( I  ,  J  - 1 )  " 1  +  SGYd.J)"1)'1  .  2 


(23) 

(24) 

(25) 


This  arrangement  occurs  because  sees  the  xx  and  xy  conductivities  in 

parallel,  while  sees  the  yx  and  yy  conductivities  in  series. 


The  conductivity  matrix  for  at  cell  (I,J)  is  analogously  described 
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a°(I-J)xx  =  (SGX(I-l.J)  *  +  SGX(I.J)  *)  ‘  .  2  (26) 

a0(I,J)  -  (SGXY(I-1 , J)  1  +  SGXY(I.J)'1)'1  •  2  (27) 

xy 


ff0(I,J)  -  (SGXY(I-l.J)  +  SGXY(I,J))/2 


(28) 


ffo(I.J)  -  (SGY(I-l.J)  +  SGY(I,J))/2  +  SGCY(I.J) 


(29) 


Note  that,  although  the  physical  conductivity  tensor  is  symmetric  at 

each  cell  (SGXY(I,J)  =  SGYX(I,J)),  the  mathematical  conductivity  just 

described  is  not  symmetric  (<70(I,J)  *  o0(I,J)  ). 

xy  yx 


The  dielectric  properties  of  the  scatterer  are  represented  by  five 
analogous  arrays,  EPX,  EPY,  EPXY,  EPCX  and  EPCY.  These  are  combined 
in  the  same  way  to  form  the  mathematical  permittivities  e^I.J)^.  at  the  E^ 


and  E 

y 


evaluation  points  of  cell 


(I.J). 


Due  to  the  anisotropic  cross-terms,  it  is  necessary  to  know  E  at  the 
E^  evaluation  points.  This  is  done  by  simple  linear  interpolation, 


EY(I,J)  -  ( EY ( I , J )  +  EY ( 1+1 , J )  +  EY ( I , J - 1 )  +  EY( 1+1 , J - 1) )/4  (30) 

E^  at  the  E^  evaluation  points,  EX(I,J)^  is  obtained  the  same  way.  The 
matrix  difference  equation  (21)  is  then  solved  twice  at  each  cell  and  each 
time  step,  once  centered  at  and  to  advance  EX(I,J),  and  once  centered  at  and 
to  advance  EY(I,J). 


The  following  notation  is  also  used: 


-[f00]'1[^o)At 

QXX(I,J)  is  the  (1,  1)  component  of  e  evaluated  at  the 

points 

-[«„]  1  [<^o  1 At 

QXY(I,J)  is  the  (1,  2)  component  of  e  evaluated  at  the  E^ 

points 


-[«„]  [o0]At 

QYX(I,J)  is  the  (2,  1)  component  of  e  evaluated  at  the  E^ 

points 


QYY(I,J)  is  the  (2,  2)  component  of  e 
points 


-  [«  ]  " 1  [  ]  At 


evaluated  at  the  E 


-[«*]’  kolAt  .I 

SXX(I.J),  SXY(I , J) ,  SYX(I.J)  and  SYY(I.J)  are  ([I]  -  e  )[a0] 

correspondingly  located  and  defined. 


Additionally,  if  we  refer  back  to  eqs .  (15)  -  (17), 

TAUXX(I.J),  TAUXY  ( I ,  J  )  ,  TAUYX(I.J)  and  TAUYY  ( I ,  J )  are  [cr0  ]  ” 1  [  -  e0] 

analogously  located  and  defined,  and 

RXX(I ,  J)  ,  RXY  ( I ,  J )  ,  RYX(I.J)  and  RYY(I,J)  are  [u0]'V0  -  ab) 

analogously  located  and  defined. 

S 

It  is  necessary  to  evaluate  both  components  of  [E  ]  as  given  in  eq. 

(11)  at  both  E  and  E  points  in  each  cell.  The  above  conventions  indicate 
\  /  x  y  r 

how  to  combine  the  ct0  and  tensors  for  an  inhomogeneous  scatterer  so  this 

complete  set  of  evaluations  may  be  achieved.  In  particular,  we  denote  E^ 

S  S  S 

as  E  evaluated  at  the  E  mesh  point  and  E  as  E  evaluated  at  the  E  mesh 
y  x  *  yx  x  y 

point . 


Thus,  ES  (I.J) 


is  , 


from  eq .  ( 11)  , 
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E®  (I.J)  -  -  TAUXX(I,J)E^C(I,J)  -  TAUXY(I,J)E*"C(I,J) 

AA  XX  XV 


-  RXX(I,J)E^C(I,J)  -  RXY(I,J)E^C(I,J) 


(31) 


where  E_  (I,J)  is  the  analytically  specified  E  field  evaluated  at  E 
XX  .  X  .  x 

mesh  points,  and  E1  (I,J)  is  the  analytically  specified  Elnc  field 

xy  S  Y 

evaluated  at  E^  mesh  points.  Additionally,  E  (I,J)  is 


Eyy  ( I  ’  J )  “  '  TAUYY(I,J)Ey"C(I,J)  -  TAUYX(I,J)E^C(I,J) 
-  RYY  ( I ,  J  )  Eyy C  ( I ,  J  )  -  RYX(I,J)E™C(I,J) 


(32) 


i.TlC  3.TIC 

where  E  (I,J)  is  the  analytically  specified  E  field  evaluated  at  the  E 

yy  •  y  v 

mesh  points  and  E  (I,J)  is  the  analytically  specified  E  field  evaluated 

yx  S  S  X 

at  the  E^  mesh  points.  Finally,  Ex  (I,J)  and  Evx(I,J)  are,  in  analogy  with 

eq.  (30), 


yx 


ES  (I,J) 
xy 

ES  (I,J) 
yx 


(Eyy(I.J)  +  E^I+l.J)  +  Eyy  ( I » J '  1 )  +  Ey^(I+l ,  J  - 1)  )/4  (33) 
(Eyv  (I  >  J)  +  Evv(I,J+l)  +  eL(I-I.J)  +  E®(I-l,J+l))/4  (34) 


It  is  also 
eq .  (9)  at  both  E 

T 

above,  we  let  J 

xy 


T 

necessary  to  evaluate  both  components  of  [J  ]  as  given  in 


and  Ey  points  in  each  cell.  Using  the  same  convention  as 

T  T  T 

be  J  evaluated  at  the  E  points  and  J  be  J  evaluated 
y  x  yx  x 
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T 

at  the  E^  points.  Then  equation  (21),  where  [J  ]  is  actually  required,  uses 
T 

[J  ]  in  the  form 


kofVh  -  -  [ES]  +  [a0]'1[Jf]  -  [ffol'V  X  Hscat] 


(35) 


S  T 

We  have  already  described  how  to  find  the  [E  ]  contribution  to  [J  ] .  It  is 

easy  to  evaluate  [o0]  because  J^.  is  a  prescribed  analytic  forcing  term 

which  can  readily  be  evaluated  at  either  the  E  or  the  E  points.  The 

scat  ^ 

troublesome  term  is  [V  x  H  ] .  This  will  have  both  an  x  and  a  y  com¬ 
ponent,  each  of  which  must  be  evaluated  at  the  E^  and  the  points. 

Let  us  designate  (V  x  HSCat)^  as  the  x-component  of  this  term 

evaluated  at  the  E  points: 

x 


(V  x  HSCat) 


XX 


HZd.Jt  -  Hzq.j-l) 
Y( J )  -  Y(J-l) 


(36) 


The  y-component  evaluated  at  an  point  is 


(V  X  HSCat) 


xy 


[(Hzq+i.j)  +  Hzq+i.j-m 

-  fHzq.j')  +  Hzq.j-iq 

2(X(I+1) 

-  X(I) ) 

(,HZq.J~)  +  HZ(I.J-l))  -  (HZd-l.J)  ±  HZd-l.J-1)) 
2 (X ( I )  -  X(I-l)) 


(37) 


The  y  component  evaluated  at  an  E  point  is 


(V  X  HSCat) 


yy 


Hzq.j)  -  Hzq-i.j) 

X(I)  -  X(I-l) 


(38) 
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and  the  x  component  evaluated  at  an  point  is 


(V  x  Hscat) 


1  f  (HZ(I .  J+l)  +  HZn-l.J+l')')  -  (HZCI.J')  +  HZfl-l.J 


2(Y(J+1)  -  Y( J ) ) 


HZa.J~>  +  HZ(I-l.Jl)  -  (HZ(IJ>1)  -  HZCI-l.J-1 


2 (Y( J )  -  Y( J - 1) ) 


Consequently,  for  example,  eq.  (21)  for  advancing  EX(I,J)  in  all  its 
glory,  becomes 


EX(I,J)n+1/2  =  QXX(I,J)EX(I,J)n'1/2  +  QXY(I,J)EY(I,J)n"1/2 


+  (1  -  QXX(I,J))E^x(I,J)n  -  QXY(I,J)E^y(I,J)n 


-  SXX(I,J)(Jf(I,J)”x  -  (V  x  HSCat)"x) 


-  SXY(I,J)(Jf(I,J)x  -  (V  x  HSCaV  ) 


where  QXX  and  QXY  are  defined  after  eq.  (30),  EY(I,J)x  is  defined  by  eq. 

(30),  E^  (I,J)  is  defined  by  eq.  (31),  E^  (I,J)  is  defined  by  eq.  (33),  SXX 
xx  xy 

and  SXY  are  defined  after  eq.  (30),  and  Jf^I,J^Xy  are  the  forcinS 

currents,  (V  X  HSCat)  is  defined  by  eq.  (36),  and  (V  x  HSCat)  is  defined 

xx  ^  —  xy 

by  eq.  (37) . 


The  scalar  equations  for  advancing  HZ(I,J),  eq.  (20),  is  much  easier  to 
implement  than  the  matrix  equation  for  advancing  (E(I,J)].  We  now  need  to 
define 

XMUZ(l.J)  as  the  bulk  permeability  at  cell  (I.J),  ^  of  eq.  (A), 

•A: 

SGMZ(I.J)  as  the  bulk  magnetic  conductivity  at  cell  (I.J)  ,  o0  of  eq. 
(A). 


QMZZ( I , J ) 


-XMUZ(I.J) 

e 


SGMZO  ,J) 


At 


(Al) 


SMZZ(I.J)  -  (1  -  QMZZ(I,J))/SGKZ(I,J)  (A2) 

TAUMZZ(I.J)  -  ^‘1(Mo0  -  Po)  (43) 

evaluated  at  the  center  of  cell  (I.J),  and 

RMZZ(I.J)  -  ^o'1(°o  -  (44) 

also  evaluated  at  the  center  of  cell  (I.J) 

The  murderously  complicated  interpolations  involved  in  advancing  E  do  not 

occur  in  advancing  H  partly  because  H  is  the  only  component  of  H  present, 

Z  2 

and  partly  because  Hz  is  evaluated  at  the  center,  not  on  an  edge  of  the 
cell . 


From  eq.  (10),  (I,J)  is  then 

HS  (I,J)  =  -  TAUMZZ(I,J)HinC(I,J)  -  RMZZ(I,J)H^C(I,J)  (A5) 

zz  ZZ  ZZ 

where  HinC(I,J)  is  the  analytically  specified  H*nC  field  evaluated  at  the 
zz  z 

cell  centers. 

★  T 

One  also  need  only  evaluate  J  (I,J)  of  eq .  (8)  at  the  cell  centers. 

*T  Z 

Equation  (20),  where  Jz  (I,J)  actually  appears,  uses  the  form 
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% 

V 


*  -1  *T 

(<*o>  J2  (I.J) 


+  (a*)'1^  x  ESCat). 


(46) 


S  * 

Here,  we  already  have  found  H  (I,J),  ^f(-*-'^)Z2  i-s  a  prescribed  magnetic 


current  density  (which  would  be  zero  on  any  physically  real  problem) ,  and 


_scat. 

(V  x  E  )  is  just 


(V  x  ESCat) 


EY(I+1,J)  -  EY ( I , J )  EX(I , J+l)  -  EX ( I , J ) 


zz  X0 (1+1)  -  X0 ( I )  ‘  Y0 (J+l)  -  Y0(J) 


(47) 


$ 


v.1 


«o  * 


Thus,  eq.  (20)  for  advancing  HZ(I,J)  becomes 

HZ (I , J)n+1  -  QMZZ(I ,J)HZ(I,J)n  +  (1  -  QMZZ(I,J))H^z(I,J)n+1/2 
-  SMZZ(I,J)(jt(I,J)n+1/2  +  (V  x  Escat  n+1/2 


Introduction  of  Frequency  Dependence 


Let  us  assume  the  frequency- dependent  term  in  Eq .  (19)  has  a  kernel  of 

the  form 


K(m) 


M 


m=l 


a  e 
=m 


■V 


(49) 


where  has  the  units  of  conductivity.  This  assumption  is  equivalent  to 
expanding  the  f requency-dependence  of  the  material's  electrical  properties 
in  a  Prony  series  under  the  constraint  that  all  the  poles  be  on  the  real 
axis.  (Appendix  1  indicates  how  one  may  relax  the  real -poles  -  only 
constraint . ) 


Equation  (19)  then  becomes 


e 

—co 


3E 


scat 


at 


+  O' 


M 

scat  .  \  .scat 

+  )  a  •  J 
[_  -m  ~m 

m=l 


(50) 


where 


jScat 

— m 


(t) 


-P  *  rt  3ESCat(t')  fit' 


J 


at' 


e  m  dt' 


(51) 


T 

In  eq.  (50),  J  is  still  given  by  eq.  (35),  but  with  the  understanding  that 

S  S 

E  has  the  frequency-dependent  term  restored.  In  other  words,  E  is  now 

represented  by  eq.  (11),  not  eqs .  (31)  and  (32). 

S  C  a  t 

The  J  of  eq .  (50)  are  not  clearly  identifiable  either  as  conduction 

TO 

or  as  displacement  currents.  We  shall  coin  the  name  "Prony  currents"  for 

scat 

them.  Equation  (50)  for  (t)  is  much  easier  to  recognize  if  we  dif¬ 

ferentiate  it  once;  its  homogeneous  solution  is  just  a  decaying  exponential: 


32 


njiT/W.v 


I 

s 


ajscat  „  aEscat 

m.  -  .  a  .scat  _z - 

at  +  -  at 


(52) 


In  the  2D  TM  case,  the  components  of  eqs .  (50)  and  (52)  then  comprise 

„„  .  ,  .  j-  -scat  ,  .scat 

2(M+1)  coupled  first-order  differential  equations  for  E  and 

Ideally,  these  equations  should  all  be  advanced  from  (n-l/2)At  to  (n+l/2)At 

simultaneously  each  cycle.  A  technique  for  doing  this  is  also  described  in 

Appendix  1 . 

However,  the  present  code  actually  implements  a  slightly  less  accurate 

algorithm  where  gscat  alone  is  advanced  first  in  each  cycle,  and  then  the 

Jscat  are  advanced  separately.  Finally,  a  correction  is  made  to  the  ad- 
— i in  . 

vanced  ESca  to  account  for  the  effects  of  the  J  . 

At  this  point,  it  is  most  instructive  to  go  back  to  eq.  (7)  and  perform 
a  rearrangement: 


V  X  HSCat  -  £ 


3E 


scat 


.scat 


3E 


inc 


■«  3 1 


+  2o  •  E  +  '  I£o)  *  at 


rt  3ET(t' ) 

+  (go  -  IO  •  ElnC  +  Jj<(t-t')  •  a-  ,  -  +  J 


(53) 


The  inhomogeneous  part  of  this  equation  can  be  written 


[jT]  +  [jP]  .  .  [o0][ES]  +  [JP]  +  [Jf]  -  [V  x  HSCat]  +  [JP] 


(54) 


where  [ E^ ]  and  [JP] 


are 


C 


Elnc] 


,  inc , 
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(55) 


n 


UP]  -  f  [K(t-t')]  [ET(t')]dt' 


with  [r]  and  [R]  still  respectively  given  by  eqs.  (15)  and  (16). 


Substituting  eqs.  (54) -(56)  back  into  (53)  then  results  in 


~  ,  scat  T  P 

ioo  *  at  +  2o  •  E  -  -  J  -  J 


This  equation  is  just  (19)  with  the  frequency-dependent  term  transferred  to 

p 

the  right  and  represented  as  J  .  It  is  advanced  according  to  eq.  (21): 


n+1/2  '^€aJ  t^oJAt  n-1/2 

[E(I , J) ]n  L/Z  =  e  [E(I,J)]n 


-  ([I]  -  e 


-[«  ]*  [^0]At 


)  [<701  1[jT(i.j)  +  JP(I,J)]n  (58) 


The  problem  with  direct  application  of  this  procedure  is  that  we  do  not 

p  •  S  O  3.  t 

know  the  portion  of  J  associated  with  E  at  nAt  until  we  have  advanced 

the  Prony  currents,  and  we  cannot,  strictly  speaking,  do  that  until  we  have 

,  ,  „scat 

advanced  E 


As  mentioned  previously,  the  code  does  not  presently  utilize  the  proce¬ 
dure  described  in  Appendix  1  for  simultaneous  advancement  of  pscat  ancj 
J^Cat.  Rather,  we  first  find  an  intermediate  value  of  Escat  obtained  with 


effects  of  the  Prony  currents  omitted: 


;E(I,J)]lnt  -  e 


«  I’Vol^t 


[  E  ( I  ,  J ) 


, n-1/2 


-  ([I] 


[£T0  ]At 

)[*o) 


1[JT(I,J)]n 


(59) 


The  procedure  for  obtaining  this  intermediate  value  is  identical  to  the 
procedure  described  in  the  previous  section  for  advancement  through  a  total 
cycle  in  the  absence  of  frequency-dependent  effects.  Its  implementation  in 
the  code  is  also  identical  to  what  was  described  in  the  previous  section. 

p 

Let  us  next  turn  to  the  advancement  of  the  total  Prony  current  J  as 
given  by  eqs .  (49)  and  (56): 


M 


m=l 


(60) 


where 


T 

sOt) 


-Bt  ft  ar<t')  Bf 


J, 


at' 


m  j  i 

e  dt ' 


(61) 


Equation  (61) ,  like  (51) ,  is  made  more  recognizable  by  differentiating  with 
respect  to  time: 


T  T 

3J  T 

_ ns  ,  a  7t  _  _ 

at  ^m- m  at 


3(ElnC  + 


ESCat) 


at 


(62) 


The  equation  for  exponential-difference  advancement  of  this  result  is 


JT(I,J)n+1/2  -  e  JT(I ,  J)n'1//2 


-  e 


Frequency  dependence  in  the  magnetic  properties  of  the  material  can  be 
treated  in  an  exactly  dual  manner  to  what  we  have  just  described  for  the 
electrical  properties. 
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At  present,  we  have  not  coded  up  magnetic  frequency  dependence,  nor 
have  we  combined  frequency  dependence  with  off-diagonal  type  anisotropy. 
Both  these  generalizations  would  be  perfectly  obvious  extensions  of  what  has 
been  done,  but  we  cannot  conceive  a  cannonical  problem  we  could  check  the 
results  against. 


Not  mixing  off-diagonal  anisotropy  and  frequency  dependence  means  we 


only  treat  diagonal  tensors.  If  eq .  (65)  is  substituted  in  eq.  (66),  u'e 

V.  t-  . 


out  a  i.  n 


-1, 


;E(I ,J) ]n+1/2  =  [E(I,J)]lnt  -  ([I]  -  e 


-l«  ]  ko]At 


M 


} 


(SAm(I,J)][J^(I,J)n+1/2  +  J^(I,J)n'1/2j 


(67) 


m=l 


where 


;sAm(i,j)]  =  [ <70 ( i , j ) ]  "(Am(i,j); 


(68) 


In  the  actual  code,  two  arrays  are  used  to  describe  the  material  ef¬ 
fects  of  each  term  in  the  Prony  series.  Let  <70(I,J)  ^  represent  the  xx 
element  of  the  inverse  of  the  matrix  described  by  eqs.  (22) -(25)  at  cell 
(I,J).  Let  ct0(I,J)  similarly  represent  the  yy  element  of  the  inverse  of 
the  matrix  described  by  eqs.  (26)-(29)  at  cell  (I,J). 

Moreover,  let  AAMX(I,J)  represent  the  xx  element  of  the  mth  Prony 
tensor  of  the  bulk  material  at  cell  (I,J),  and  let  AAMY(I,J)  represent  the 
corresponding  yy  element.  Then  the  xx  element  of  [SA  (I,J)]  which  actually 
relates  the  x  component  of  JT(I,J)  to  the  x  component  of  E(I,J)  is  called 


(69) 


SAMX(I.J)  =  a0  (I , J)  ~  (AAMX(I.J-l)  +  AAMX(I,J))/2 


Similarly,  the  yy  element  of  [SA  (I,J)]  which  actually  relates  the  y  com¬ 
ponent  to  the  y  component  of  E(I,J)  is  called 


SAMY(I.J)  =  o0(I,J)"^  (AAMY(I-l.J)  +  AAMY(I,J))/2 


(70) 


In  keeping  with  our  simplification  of  not  mixing  off-diagonal  anisotropy 
with  frequency  dependence,  we  ignore  any  possible  off-diagonal  nonzero 
values  in  [ SA  (I,J)1. 

It  turns  out  that  only  SAMX(I,J)  and  SAMY(I,J)  need  actually  be  stored. 

That  is,  it  is  not  necessary  to  assign  arrays  for  keeping  a0 (I , J) ' ^ , 

-1  xx 
£7  o  ( I ,  J  )yy  ’  ^A>tX(I,J)  and  AAMY(I ,  J )  . 

Consequently,  the  actual  equation  used  in  the  code  for  implementing  the 
x-component  of  the  Prony  correction  is 


EX(I,J)n+1/2  =  EX(I , J) lnt 


M 

-  (1  -  QXX ( I , J ) )  ^  SAMX(I , J)  XJMSX(I,J)n  (71) 

m=l 


where 


XJMX( I , J ) n  =  [JT(I , J)n+1/2  +  JT(I ,J)n'1/2]/2 


(72) 


Similarly,  the  actual  equation  used  for  implementing  the  y-component 
Prony  correction  is 


EY(I,J)n+1/2  -  EY(I,J)lnt 


M 

-  (1  -  QYY(I.J))  ^  SAMY(I.J)  XJMY(I,J)n 
m=l 


where 


XJMY(I , J)n  -  [JT(I,J)n+1/2  +  JT(I,J)n'1/2l/2 


of  the 


(73) 


(74) 


Transformation  from  the  Near-field  Time  Domain  to  the 
Far- field  Frequency  Domain 


The  foregoing  work  described  the  determination  of  the  scatterer’s 
electromagnetic  response  and  associated  near  fields.  We  are  actually  inter¬ 
ested  in  the  RCS ,  which  is  a  far-field  quantity.  Now  we  shall  describe  how 
the  code  extracts  the  RCS  from  the  near- field  results.  In  this  process,  we 
also  transform  from  time  domain  to  frequency  domain. 


Any  electromagnetic  field  can  be  expressed  in  terms  of  an  electric  and 
magnetic  vector  potential,  A  and  A  .  These  vector  potentials  (in  the 


frequency  domain)  obey  the  inhomogeneous  wave  equations. 


V2A  +  k2A  «  -p0J 


(75) 


„  *  „  *  * 
V2A  +  k2A  «  e0J 


(76) 


Here  J  is  the  fictitious  magnetic  current  density  often  found  useful  in 
manipulating  Maxwell's  equations.  Equations  (75)  and  (76)  can  be  general¬ 
ized  to  apply  to  any  linear  medium,  although  we  shall  find  their  free-space 
form  adequate  for  our  uses . 


In  2D,  we  define  the  far  field  to  be  the  region  where  all  fields  drop 
-1/2 

off  as  r  '  ;  i.e.,  where  all  the  faster  falling  terms  have  vanished.  We 

can  then  separate  the  electric  and  magnetic  fields  into  two  parts, 


E  +  E 
-e  ~m 


(77) 


H  +  H 
— e  — m 


(78) 


AO 


.V.  -• 


M 


1  \.' 


- 

S  , 


v.  / 


■i 


The  e - subscripted  parts  are  associated  with  the  electric  vector  potential  A, 
and  the  m- subscripted  parts  are  associated  with  the  magnetic  vector  poten- 

"$C 

tial  A  .  In  particular,  if  we  call  Vq  and  Z0  the  admittance  and  impedance 
of  free  space,  we  can  show  in  the  far  field  that 


He  =  V  x  A /nQ  «  iwY0ir  x  A 


-k  ~k 

E  =  V  x  A  / f n  =  i^Z0i  X  A 
- m  -  /  o  o-r  - 


iwA^  =  iw(i,A,  +  i  A  ) 
~t  <f>  <f>  — z  z' 


“k  "k 

-  iw(i  , A ,  +  i  A  ) 
<f>  <p  ~  z  z 


Here,  a  t  subscript  indicates  that  only  the  transverse  components  ($  and  z) 
are  retained.  Equations  (79) -(82)  are  analogous  to  3D  formulas,  and  depend 
on  the  fact  that 


i (kr-wt) 


is  a  valid  far- field  frequency- domain  2D  solution  of  the  wave  equation  even 
if  the  more  general 

f  ( t _ - _ r/c)  /  o  /  \ 


is  not  a  valid  time-domain  solution.  These  equations  tell  us  that  if  we  can 
evaluate  A  and  A  ,  we  can  find  the  2D  RCS  without  undue  complication. 


In  two  dimensions,  the  Green's  function  for  the  scalar  wave  equation 
in  the  frequency  domain  obeys 


V2G(r|r')  +  k2G(r|r')  =  S(r-r') 


(85) 


where  r  is  the  scatterer  location  and  r'  is  the  observer  location.  This 
equation  has  solution 


G(r|r')  -  -  ^H^CkR) 


(86) 


Thus,  at  least  in  cartesian  coordinates,  A  and  A  become 


A(r' ,«) 


rfiH(1)(kR) 

IJ-fc— 


J(r ,w)dr 


(87) 


A*(r ' ,w)  - 


rriH(1) (kR) 


1  (X.«)d£ 


(88) 


For  the  far  field  region,  G(r|r')  asymptotically  approaches 


G(r  |  r '  ) 


(89) 


This  expression  may  be  further  manipulated  by  letting  r'  replace  R  in  the 
denominator  of  the  radical.  The  phase  term  requires  a  bit  more  care: 


kR  -+  k(r '  -  i' 


r)  —  kr'  -  (kx  cos<f>'  +  ky  sin^') 


(90) 


where 


ir,  =  J^COS^'  +  iy  sin<£'  (91) 

is  a  unit  vector  pointing  from  the  target  to  a  far-field  observer. 

Using  these  expansions,  we  can  rewrite  the  formula  for  A  in  the  far 
field  as 


A(r' ,co)  = 


Moe 


-3i?r/4 


/  2  ikr'  ffT,  .  -ik(x  cos6'  +  v  sind)  ,  .... 

v  jrkr7  e  JJ-(£,w)e  d£  (92) 


A  corresponding  expression  exists  for  A  (r' ,u) .  The  far-field  expression 


for  E  then  becomes 
— e 


Ie(r',w)  =  iw  At(r' ,w)  = 


Wfi0e 


■  3i7r/4 


7  elkr'  U  J,(r,u,)e-ik<*  +  S'  sl"*'>dr 


(93) 


Similarly,  H  becomes 


H  (r' ,w)  -  iwA  (r' ,«)  = 


4 


(94) 


xco  c  ne 


wkr' 


ikr 


'  JJ  J*( r,o,)e-ik(x  cos^'  +  y  sin*'>dr 


Analogous  formulas  exist  for  E  and  H  . 
°  — m  “e 


Equations  (93)  and  (94)  are  not  directly  applicable  to  the  output  of 
our  2D  Maxwell  solver  as  these  equations  demand  the  frequency-domain  J  and 
J  ,  while  the  Maxwell  solver  outputs  the  time-domain  currents. 


Let  us  say  we  want  E  (r*  w  )  where  w 
J  — e  —  ’  q'  q 

of  interest.  We  can  then  write 


is  one  of  N  discrete 

q 


frequencies 


E  (r' ,w  )  - 
-e  ~  '  q' 


- 3in/U 


(95) 


where  k'  is  the  wavenumber  pointing  towards  the  observer  at  to  and  where  we 
q  q 

have  interchanged  the  order  of  time  and  space  integration  after  replacing 

Jt(r,w)  with  its  inverse  Fourier  representation. 


Analogously , 


H  ( r a)  )  becomes 
— m  q 


lco  e „e 


-  3i7r/4 


H  (r'  ,w  ) 

— m  x  —  ’  *-%  ' 


ft  ir 


\  > 


ik  r' 

,  q 


J^e  -  dt  JJ 


Equations  (95)  and  (96) ,  and  the  two  companion  equations  for  and  Hg 
are  in  a  form  which  is  compatible  with  our  time-domain  Maxwell  solver. 
Taking  the  time  integration  outside  the  space  integration  is  vitally  impor¬ 
tant  to  the  efficiency  of  our  algorithm.  Were  this  not  done,  J.t(r,t)  and 
J  (r,t)  would  have  to  be  Fourier  transformed  at  every  point  before  being 
integrated  over  space.  The  form  of  eqs .  (95)  and  (96)  replaces  this  enor¬ 
mous  computation  with  a  single  Fourier  transform  on  the  result  of  the  space 
integral . 

If  we  let  JSCat  represent  the  total  current  (conduction,  displacement 
and  Prony)  associated  with  the  scattered  electromagnetic  field  (see  eq. 
(53)) 


Tscat 
J  -  e 


«  3t 


„scat  ,  ,  T  . 

+  2o  *  I  +  <«  '  Ie0) 


3ET ( t ' ) 


+  <£o  -  ICTb)  •  1  +  JroS<t-t:')  *  — ^ ; 


and  if  we  substitute  J  for  in  eq.  (95),  of  eq.  (95)  becomes  the 

scattered  field  ESCat  of  the  first  section  unless  magnetic  materials  are 

_ .  rr.1 .  '  TSCat  "  _ . ,  _  .  1 .  . _  _ _ _ .* , 


present.  That  is,  J“ 


integrated  over  the  scatterer  cross-section  accord¬ 


ing  to  eq .  (95)  gives  gscat  in  the  absence  of  magnetic  materials. 


RCS(w  )  —  2n 

q 


-scat ,  ,  .  /  ,  2 

E  (r',w  )  •  Jr' 

—2 - y - 

Ei(r ' ,«  ) 


In  this  convention,  A  and  E  are  zero. 


However,  it  is  possible  to  replace  the  area  integral  of  eqs.  (95)  and 
(96)  with  a  contour  integral  by  means  of  Huygens  principle.  In  particular, 
let  S  be  a  closed  contour  which  completely  surrounds  the  target.  (For 
instance,  let  S  be  a  rectangle  defined  by  x  -  xx ,x2  and  y  -  y1(y2.)  Let 

S  C S  u  SC s t 

E  and  H  be  the  scattered  fields  which  our  Maxwell  solver  predicts 
will  exist  on  S  due  to  the  time-domain  illumination.  Let  n  be  an  outward¬ 
pointing  unit  normal  on  S .  If  we  then  remove  the  scatterers  and  its 
currents,  but  let  an  electric  surface  current 


K  -  n  x  H 


scat 


(99) 


and  a  magnetic  surface  current 


K  *=  -  n  x  E 


scat 


(100) 


flow  on  S,  the  scattered  electromagnetic  field  will  be  replicated  outside  S. 
This  means  the  area  integral  of  eq.  (95)  may  be  replaced  by 


S  (n  x  HSCat(r  ,t))  e  ^  ^  As  -  ISCat(t), 

p  ~ p  "t  p q  v  '\ 


(101) 


where  the  summation  over  p  represents  integration  over  the  finite  difference 
cell  edges  which  lie  on  S.  This  summation  is  represented  by  I  because  it 
has  the  dimension  of  amperes.  In  the  actual  code,  it  goes  by  the  name 
XIEST  . 


Analogously,  the  area  integral  present  in  eq.  (96)  for  H^(r  .t^)  can  be 


£  (-n  x  ESCat(r  ,t))  e 
p  —  ~p  t 


i 


(102) 


(Remember  the  t  subscript  on  I,  X*,  or  anything  else  implies  evaluation  with 
the  r  component  omitted.)  In  the  code,  this  variable  is  called  XIMST^. 

scat 

Equations  (95)  and  (101)  indicate  the  z-component  of  H  is  as¬ 

sociated  with  the  ^-component  of  E  (r' ).  The  ^-component  of  ESCat  in  eq. 

e  q 

(102)  has  the  same  association,  as  the  cross  with  jL  in  eq.  (80)  proves. 
These  scattered  fields  are  the  TM  solution. 


scat 

Analogously,  the  ^-component  of  H  in  eq.  (101)  and  the  z-component 

scat 

of  E  in  eq.  (102)  relate  to  the  (decoupled)  TE  solution,  which  we  are 

not  treating  in  detail  at  this  time. 


scat  in 

Let  us  use  X  (t)  to  denote  this  "current"  evaluated  from  the 

^  scat  ^n"t"l/2 

finite-difference  code  at  nAt.  Let  us  analogously  denote  X  (t)^  '  * 

As  the  finite-difference  calculation  progresses,  we  can  then  keep  running 

summations  for  each  frequency  w^, 


„  n  _  iw  mAt 

QscaC(t)"  -  2  I®  (t)“  e  q  At 

‘l  c  H  c 


(103) 


.  ,  /0  n  _  ,  ,,  jw  (m+l/2)At 

2scat(t)*n+1/2  _  s  Iscat(t)*m+l/2  eJ  q 

9  m«1  9 


(104) 


When  the  time-domain  finite  difference  calculation  is  complete,  these  Q's 
will  then  respectively  contain  quantities  which  are  directly  proportional  to 
the  electric  and  magnetic  contribution  to  the  RCS  at  u>  .  Note  that  it  is 
only  necessary  to  backstore  2N^  complex  quantities  during  the  time-domain 
finite -difference  calculation  in  order  to  preserve  all  the  information 
necessary  to  generate  a  monostatic  RCS  as  a  function  of  w. 


*  *  *  •  ’  /  *  H  «  '  j  ' 


47 


The  symbol  Q  is  used  in  eq.  (103)  because  it  represents  a  quantity  with 

units  of  coulombs.  In  the  code,  it  is  written  SIEST  .  The  scattered 

scat  ^ 

electric  field  associated  with  Q  (tl  is 

q  t 


pScat .  ,  v 

E  (r' ,u  ) 
_e  q 


8tt 


-  3i7r/4 


ik  r' 

q  Ascat, 

;  Q  (“) 
^q  v  ' 1 


(105) 


The  code  tracks  the  ^-component  of  this  quantity  as  EPES^,  when  J r'  is 
normalized  out.  Analogously,  the  scattered  magnetic  field  associated  with 

Qq  (t)t  IS 


-  3i?r/4 


..SC3t  /  .  x 

H  (r' ,w  ) 
— m  q 


rr  ik  r' 

/  2  a  .scat  .  . *» 

;  r  e  M  Q  (“>) 

?rk  r '  '  '  t 

q 


(106) 


In  the  code,  the  z- component  of  this  quantity  is  HZMS^  when  Jr’  is  normal¬ 
ized  out.  The  scattered  electric  field  associated  with  the  magnetic  current 
is  obtained  by  combining  eqs .  (80)  and  (106), 


-scat.  ,  . 

E  (r* ,w  )  - 
~m  q 


-3i?r/4 
lw  e  ' 


f~Z  ik  r'  ,  .  , 

~T~~7  e  «  i  x  Qscat(<o) 
ffk  r'  r  q  7 

q 


(107) 


In  the  code,  the  ^-component  of  this  quantity,  less  the  Jr’  factor,  is 


EPMS  . 

q 


The  radar  cross  section  then  becomes 


RCS (w  )  -  In 

q 


(fat(r,  ,»q)  +  E^cat(r',^o))yr' 

E1  (r '  ,u>_) 


(108) 


In  the  code,  SCATXE  is  the  EScat  contribution  to  the  ratio  inside  the 

q  e  ’  ' 

scat 

absolute  value  signs,  and  SCATXM  is  the  E  contributions. 

q  m 


APPENDIX  1.  SOLUTION  OF  THE  TENSOR  WAVE  EQUATION 


Let  us  consider  the  case  where  anisotropy,  but  not  frequency  depend 
ence,  is  present, 


[<J  DIE]  +  [a0][E]  =  [J] 


(A1 


This  matrix  differential  equation  has  a  homogeneous  solution 


[E] 


h 


-uj'VcJt 

e  [A] 


(A2 


and  a  particular  solution 


[  E )  p  =  [  cr0  ] " 1  [  J  ] 


(A3 


giving  a  general  solution 


[E] 


-lej'VoJt  .j_ 

e  [A]  +  [a0 ]  i[J) 


(M 


The  constant  vector  [A]  may  be  evaluated  at  (n  -  l/2)At: 


[  E  ] n" 1/2  =  [A]  +  [<70  ]  ~ 1  [  J  ] n 


(AS 


This  gives  the  new  E-field  vector  in  terms  of  the  old, 


